icr_mcmc::icr_mcmc(long inN,  long in_niter, string infilename, string outfilename){
   //constructor
   rands = new icr_rands;
   N = inN;
//cout<<N<<endl;
   niter = in_niter;
//cout<<niter<<endl;
   filename = infilename;
   outputname = outfilename;
   parse_csv();
   set_priors();
   set_initial_values();
   //initialize output file
   //outfile = new ofstream(outputname.data(), ofstream::binary);
   //fout = new ofstream(outputname.data());
   fout.open(outputname.data(), ofstream::binary);
   double   buffer = niter;
   fout.write((char *)(&(buffer)),sizeof(buffer));
   buffer = N;
   fout.write((char *)(&(buffer)),sizeof(buffer));
   buffer = ncol;
   fout.write((char *)(&(buffer)),sizeof(buffer));

}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::run(void){
  cout<<"ICR_MCMC::RUN running "<< niter <<" iterations..."<<endl;
  /*double buffer = 6;
  for(int i=0;i<10;i++){
  fout.write((char *)(&(buffer)),sizeof(buffer));
  }*/


  for(int i = 0;i<niter;i++){
      update_mustar_sigstar();
      stickbreaking_updates();
      update_indicators();
      block_find();
      update_params();
      update_alpha();
      if(state->iternum % 10 == 0 && state->iternum >0.3*niter){write_output();}
      state->iternum++;
      //cout<<"In debugging mode so not doing all iterations"<<endl;
      //exit(1);
  }
  time_t rawtime;
  struct tm * timeinfo;
  time (&rawtime);
  timeinfo=localtime(&rawtime);
  cout<<"ICR: Finish at time: "<<asctime(timeinfo)<<endl;

  
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::find_nlines(void){
 ndata = 0;
 string line;
 ifstream myfile(filename);
 if(myfile.is_open()){
    while(!myfile.eof()){
       getline(myfile, line);
       ndata++;
    }
 } else{
   cout<<"couldn't open file: "<<filename<<endl;
 }
 myfile.close();

}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::find_ncol(void){

 ncol = 0;
 ifstream myfile(filename);
  string                line;
    getline(myfile,line);

    stringstream          lineStream(line);
    string                cell;
    rowvec tmprowvec;
    while(getline(lineStream,cell,','))
    {       ncol++;    }
 myfile.close();
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::parse_csv(void){
  find_nlines();  
  find_ncol();  
  cout<<"Number of datapoints: "<<ndata<<endl;
  cout<<"Number of dimensions: "<<ncol<<endl;
  data.zeros(ndata, ncol);
  rowctr = 0;
  ifstream myfile(filename);
  if(myfile.is_open()){
    while(!myfile.eof()){
       getNextLineAndSplitIntoTokens(myfile);
       rowctr ++;
    }

  }
  myfile.close();
//  data.print("data:");

}

////////////////////////////////////////////////////////////////////////
void icr_mcmc::getNextLineAndSplitIntoTokens(istream& str)
{
    //vector<string>   result;
    string                line;
    getline(str,line);

    stringstream          lineStream(line);
    string                cell;
    rowvec tmprowvec;
    long ctr=0;
    while(getline(lineStream,cell,','))
    {
        data(rowctr,ctr)  = atof(cell.c_str());
        ctr++;
    }
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::set_priors(void){
  L = ncol;
  h =  median(data, 0);
  rowvec maxes(ncol);
  rowvec mins(ncol); 
  for(int i =0; i<ncol; i++){
      maxes[i] = max(data.col(i));
      mins[i] = min(data.col(i));
  }
  r = maxes - mins;
//  maxes.print("max:");
//  mins.print("min:");
  //r.print("r:");
  H.zeros(ncol, ncol);
  for(int i=0;i<ncol;i++){
     H(i,i) = r[i]*r[i]/16;   
  }
  //H.print("H:");
  //priors for m
  a_m = h;
  B_m = H;
  B_m_inv = pinv(B_m);
  //priors for V
   a_V = 2*(L + 2);
   B_V = (a_V - L - 2)*H;
  //priors for S
   a_S = 2*(L + 2);
   nu = 2*(L+2); 
   B_S = (nu - L -2)*H/a_S;
   B_S_inv = pinv(B_S);
  //prior for alpha
   a_a = 10;
   a_b = 1;
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::set_initial_values(void){
   state = new icr_mcmc_state;
   state->iternum = 0;
   state->nu = nu;
   state->alpha = a_a/a_b;
   state->m = h;
   //state->m.print("m:");
   //h.print("h:");
   state->V = B_V;
   state->V_inv = pinv(state->V);
   state->S = B_S;
   //set p with stickbreakings
   stickbreaking_first();

   //set the intiial indicator variables
   state->l.set_size(ndata);
   rowvec cumsum0 = cumsum(state->p);
   cumsum0 /= sum(state->p);
   for (int i =0;i<ndata;i++){state->l(i) = rands->random_index(cumsum0); }

   state->l_rev.resize(N);
   state->njstar.set_size(N);
   //state->l.print("l:");

   state->mu_star.set_size(N, ncol);
   rmvnorm mustar_gen= rands->mvnorm_gen( state->m , state->V);
   for(int i = 0;i<N;i++){
     state->mu_star.row(i) = mustar_gen.draw();
   }
//   state->mu_star.print("mu_star:");

   state->sigma_star.set_size(ncol, ncol, N);
   riwish sigmastar_gen= rands->riwish_gen( state->nu , state->S);
   for(int i = 0;i<N;i++){
   state->sigma_star.subcube(span::all,span::all,span(i,i)) =sigmastar_gen.draw();
//     sigmastar_gen.draw().print("submat");;
     //state->sigma_star(span(),span(),span(i,i)).print("submatr");
   }
    // state->sigma_star(span(),span(),span(0,10)).print("submatr");
   //median(state->mu_star, 0).print("median mu_star:");
   block_find();
}

////////////////////////////////////////////////////////////////////////
void icr_mcmc::block_find(void){
   //fill in vector of all data point indices that correspond
   // to a given component
   for(int i=0;i<N;i++){state->l_rev[i].clear();}
   state->njstar.fill(0);
   for(int i=0;i<ndata;i++){
       state->l_rev[state->l(i)].push_back(i);
   }
   //count up number of non-zero components and populate njstar 
   // which is number in each component
   state->nstar = 0;
   long buffer;
   for(int i=0;i<N;i++){
       buffer = (long) state->l_rev[i].size();
       state->njstar(i) = buffer;
       if(buffer != 0) state->nstar ++;
   }

 //  state->njstar.print("njstar");
/*cout<<"size:\t"<<state->l_rev[5].size()<<endl;
   for(int i=0;i<(int)state->l_rev[5].size();i++){
       cout<<state->l_rev[5][i]<<endl;
   }*/
//cout<<state->nstar<<endl;
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::stickbreaking_first(void){
   state->p.zeros(N);
   double randgam0, randgam1, randbet, prod=1;
   logpn = 0;
   gamma_generator part0= rands->rgamma_gen(1);
   gamma_generator part1= rands->rgamma_gen(state->alpha);
   for(int i=0; i<(int) (state->p.n_elem-1); i++){
      randgam0 = part0();//the alpha parameter part
      randgam1 = part1();//the beta parameter part
      randbet = randgam0/(randgam0+randgam1);
      state->p[i] = prod*randbet;
      prod *= (1-randbet);
      logpn += log(1-randbet);
   }
   if(logpn < -160) logpn = -160.;
   state->p[state->p.n_elem -1] = exp(logpn);
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::stickbreaking_updates(void){
   state->p.zeros(N);
   double randgam0, randgam1, randbet, prod=1;
   logpn = 0;
   double cumsum=0;
   for(int i=0; i<(int) (state->p.n_elem-1); i++){
      cumsum += state->njstar(i);
//cout<<cumsum<<"\t"<<ndata<<endl;
//cout<<1+state->njstar(i)<<endl;
      randgam0 = rands->rgamma_gen(1 + state->njstar(i))();//the alpha parameter part
//cout<<state->alpha + ndata - cumsum<<endl;
      randgam1 = rands->rgamma_gen(state->alpha + ndata - cumsum)();//the beta parameter part
      randbet = randgam0/(randgam0+randgam1);
      state->p[i] = prod*randbet;
      prod *= (1-randbet);
      logpn += log(1-randbet);
   }
   if(logpn < -160) logpn = -160.;
   state->p[state->p.n_elem -1] = exp(logpn);
}


////////////////////////////////////////////////////////////////////////
void icr_mcmc::update_indicators(void){ //update the l[j]'s
    mat likelimat; likelimat.set_size(ndata, N);
    for(int i=0;i<N;i++){
        likelimat.col(i) = dmvnorm(data, state->mu_star.row(i),
                     state->sigma_star.subcube(span::all,span::all,span(i,i)));
 //       likelimat.col(i).print("likelimat");
//cout<<"i\t"<<i<<"\t "<<state->iternum<<endl;
 //    state->sigma_star.subcube(span::all,span::all,span(i,i)).print("sigmajstar");
    }
    for(int i=0;i<(int)likelimat.n_elem;i++){
        //check for really bad roundoff
        if(likelimat(i) < 1e-99) likelimat(i)=1e-99;
        if(likelimat(i) != likelimat(i)) likelimat(i)=1e-99;
    }
    rowvec likeli, probs, cumsum0;
    for(int i=0;i<ndata;i++){
        likeli = likelimat.row(i) - max(likelimat.row(i));
//likeli.print("likeli:");
        probs = state->p % exp(likeli); //element-wise multiplication

   cumsum0 = cumsum(probs);
   cumsum0 /= sum(probs);
//cout<<sum(probs)<<endl;
//if(sum(probs) != sum(probs)){
 //  likelimat.print("likelimat:");
//}
//state->p.print("p:");
//cout<<"sum of p"<<sum(state->p)<<endl;
//cumsum(state->p).print("cumsum p:");
//likelimat.print("likelimat:");
   state->l(i) = rands->random_index(cumsum0); 
   }

//probs.print("probs:");
//cout<<max(likeli)<<endl;
//cout<<max(probs)<<endl;
//exit(1);
    
      
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::update_alpha(void){
    double beta = a_b - logpn;
    state->alpha = rands->rgamma_gen(N+a_a-1)()/beta;
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::update_params(void){
  // update m

  rowvec sum0;
  sum0.set_size(ncol);
  sum0.fill(0.0);
 for(int i=0;i<N;i++){      
      if(state->njstar(i) !=0){
          sum0 += state->mu_star.row(i);
      }
  }

  rowvec mustar_tilde = sum0/state->nstar;

  rowvec a_m_new =  (a_m*B_m_inv + state->nstar*mustar_tilde*state->V_inv)*
                         pinv(B_m_inv + state->nstar*state->V_inv );

  mat B_m_new = pinv(B_m_inv + state->nstar*state->V_inv);
  state->m = rands->mvnorm_gen( a_m_new , B_m_new).draw();

  // update V
  double new_nu;
  mat new_S;
  new_nu = a_V +state->nstar;
  new_S = B_V;
  double thisnjstar;
  for(int i=0;i<N;i++){
     thisnjstar = state->njstar(i);
     if(thisnjstar !=0){
         new_S += trans(state->mu_star.row(i) - state->m)*(state->mu_star.row(i) - state->m);
     }
  }
  state->V = rands->riwish_gen( new_nu , new_S).draw();
  state->V_inv = pinv(state->V);
 
  // update S
if(state->iternum ==0){
cout<<"Be sure that all matrices in S and V updates are the proper inverse/non-inverse state, make sure all multiplication orders are ok."<<endl;}
  double new_nu2;
  mat new_BS;
  new_nu2 = a_S + nu*state->nstar; 
  new_BS = B_S_inv;
  for(int i=0;i<N;i++){  //the checking here is replicated and unneeded, but does it really cause a slow down?
     if(state->njstar(i) !=0){
         new_BS += state->sigma_star.subcube(span::all,span::all,span(i,i));
     }
  }
  new_BS = pinv(new_BS);
  //needs an inverse because riwish expects inverse of scale matrix...
  state->S = rands->riwish_gen( new_nu2 , inv(new_BS)).draw_noninv();
 
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::update_mustar_sigstar(void){
   rowvec new_mu_mean;
   mat new_mu_sigma,this_mat,z_tilde, sigmajstar_inverse, new_S;
   uvec these_rows;
   mat tmp;
   for(int i=0;i<N;i++){
      if(state->njstar(i) == 0){
//cout<<i<<endl;
//state->m.print("state->m");
//state->V.print("state->V");
//state->mu_star.row(i).print("mu_star:");
        state->mu_star.row(i) = rands->mvnorm_gen( state->m , state->V).draw();
//state->mu_star.row(i).print("mu_star:");
//exit(1);
        state->sigma_star.subcube(span::all,span::all,span(i,i)) =
           rands->riwish_gen( state->nu , state->S).draw();
      } else { 
          these_rows.clear();
          these_rows.set_size(state->l_rev[i].size());
          for(int j=0;j<(int)state->l_rev[i].size();j++){
             these_rows[j] = state->l_rev[i][j];
          }
//these_rows.print("these_rows:");
          this_mat =  data.rows(these_rows);
//cout<<"checking"<<endl;
//for(int ij=0;ij<state->l_rev[i].size();ij++){
//    data.row(these_rows[ij]).print("data.rows(ij):");
//    //cout<<state->l[these_rows[ij]]<<endl;
//}
//this_mat.print("this_mat:");
//exit(1);

          //this_mat.print("this_mat");
          z_tilde = mean(this_mat, 0);
          //z_tilde.print("z_tilde:");

          sigmajstar_inverse = (state->sigma_star.subcube(span::all,span::all,span(i,i)));
//sigmajstar_inverse.print("sigmajstar_before");
          sigmajstar_inverse = pinv(sigmajstar_inverse);

          //update mu_j *
          new_mu_mean =
            (state->m*state->V_inv + state->njstar[i]*z_tilde*sigmajstar_inverse)*
                 pinv( (state->V_inv + state->njstar[i]*sigmajstar_inverse));
//          new_mu_mean.print("new_mu_mean");
          new_mu_sigma = pinv(state->V_inv + state->njstar[i]*sigmajstar_inverse);
          state->mu_star.row(i) = rands->mvnorm_gen( new_mu_mean , new_mu_sigma).draw();
          //update sigma_j*
          new_S = state->S;

          for(int j=0;j<state->njstar[i];j++){
              tmp = trans(state->mu_star.row(i) - this_mat.row(j))*
                      (state->mu_star.row(i) - this_mat.row(j));
           //tmp.print("tmp");
              new_S += tmp ;
          }
 //         cout<<state->njstar[i]<<"\t"<<state->l_rev[i].size()<<endl;
          state->sigma_star.subcube(span::all,span::all,span(i,i)) =
           rands->riwish_gen( state->nu + state->njstar[i] , new_S).draw();
//cout<<state->nu +state->njstar[i]<<endl;
//(state->sigma_star.subcube(span::all,span::all,span(i,i))).print("sigmajsar_after:");
//new_S.print("new_S");
//state->S.print("state->S:");
//if( state->iternum > 3) exit(1);
//pinv(new_S).print("inv new_S");
//exit(1);
      }
   }
}
////////////////////////////////////////////////////////////////////////
void icr_mcmc::write_output(void){
//cout<<state->iternum<<endl;
   double buffer;
   buffer = state->iternum;
   fout.write((char *)(&(buffer)),sizeof(buffer));
    for(int i=0;i<N;i++){
       for(int j=0;j<ncol;j++){
          buffer = state->mu_star(i,j);
          fout.write((char *)(&(buffer)),sizeof(buffer));
       }
    }
    for(int i=0;i<N;i++){
//          mat tmp = state->sigma_star.subcube(span::all,span::all,span(i,i));
//          tmp.print("tmp");
       for(int j=0;j<ncol;j++){
          for(int k=0;k<ncol;k++){
          buffer = state->sigma_star(j,k,i);
          fout.write((char *)(&(buffer)),sizeof(buffer));
          }
       }
    }
    for(int i=0;i<N;i++){
          buffer = state->njstar[i];
          fout.write((char *)(&(buffer)),sizeof(buffer));
    }
          buffer = state->alpha;
          fout.write((char *)(&(buffer)),sizeof(buffer));
}
////////////////////////////////////////////////////////////////////////
icr_mcmc::~icr_mcmc(){
  //destructor
}
